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We study the dynamics of a Bose-Einstein condensate in a one-dimensional optical lattice in the 
limit of weak atom-atom interactions, including an approximate model for quantum fluctuations. 
A pulsating dynamical instability in which atoms periodically collect together and subsequently 
disperse back into the initial homogeneous state may occur in the time evolution. We take into 
account the quantum fluctuations within the truncated Wigner approximation. We observe that 
the quasiperiodic behavior still persists for a single realization that represents a typical experimental 
outcome, but ensemble averages show various manifestations of quantum fluctuations. Quantum 
effects become more prominent when the effective interaction strength is increased. 
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I. INTRODUCTION 

The superfluidity of a Bosc-Einstcin condensate (BEC) 
in an optical lattice has been drawing considerable atten- 
tion in the last several years Q. As is well known, super- 
flow of the BEC in free space suffers from an instability 
when the flow velocity reaches a critical value. Such an 
instability, known as Landau or energetic instability, ex- 
ists when the superfluid flow is not a local minimum of 
energy and the system may lower its energy by emitting 
phonons 0. In an optical lattice, in addition to energetic 
instability, the BEC may also exhibit dynamical or mod- 
ulational instabilities that have been the subject of much 
experimental and theoretical research over recent years 
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When the system is in the dynamically unstable regime, 
small perturbations grow exponentially in time resulting 
in irregular dynamics, loss of coherence, or an abrupt 
stop of the transport of the atom cloud 0, [1] . 

In this paper we study dynamical instabilities of atoms 
in an optical lattice for the case of weak atom-atom in- 
teractions and also taking into account quantum fluc- 
tuations of the atoms. We recently reported that, 
by appropriately selecting the strength of the atom-atom 
interactions, the corresponding classical system may ex- 
hibit a pulsating dynamical instability in which the atoms 
nearly periodically collect to a peak in lattice occupation 
numbers, and subsequently disperse back to (very close 
to) the initial unstable state. This is different from the 
conventional view, valid at strong interatomic interac- 
tions, that dynamical instabilities for BECs in optical 
lattices are associated with irregular dynamics. The dy- 
namics of an integrable double-well system provides a 
qualitative explanation of the pulsating phenomenon [20| : 
Although the instability is a result of the interplay be- 
tween the lattice discreteness and the nonlinearity that 
makes the lattice non-integrable, the dynamics of the lat- 
tice with many sites approximates the dynamics of an 
integrable system. Related classical pulsations starting 



from already compressed atom distribution in a lattice 
have been discussed in within the framework of the 
nonpolynomial Schrodinger equation. 

As the emergence of dynamical instabilities can also 
be closely related to various fundamental effects in 
noncquilibrium quantum dynamics, for example dissipa- 
tivc transport [l5[ and quantum phase transitions [2ll |. 
it is particularly interesting to address how an entirely 
classical description of the pulsating instability [l^l , em- 
bodied in the Discrete Nonlinear Schrodinger equation 
(DNLSE), is modifled as a result of quantum fluctua- 
tions. In this paper we provide an approximate analysis 
of the quantum effects on the pulsating instability, as 
well as amend the discussion of Ref. [1^ with additional 
details and angles. 

We consider different dynamically stable initial sta- 
tionary superfluid flow states that are instantaneously 
transferred to a dynamically unstable regime by chang- 
ing the strength of atom-atom interactions. We incorpo- 
rate quantum fluctuations of the atoms using a stochas- 
tic phase-space method, the Truncated Wigner Approx- 
imation (TWA), [13, M, H m, i^. We flnd that the 
quasiperiodic behavior is still observable in individual 
stochastic realizations that represent typical outcomes 
of individual experiments. In fact, the damping of the 
quasiperiodic pulsating instability in a single stochastic 
trajectory remains very weak even in the presence of en- 
hanced quantum fluctuations and for the case of an ini- 
tial state with a substantial noncondensate atom fraction. 
However, as the timing, shape and location of the pul- 
sation events in each realization change due to quantum 
effects, in the quantum mechanical ensemble average the 
wave function revival become progressively weaker, with 
an increased damping rate, when the effective interaction 
strength is increased. Other ensemble averages provide 
additional quantitative information about the effects of 
quantum fluctuations. For instance, in the limit of weak 
quantum fluctuations, the amplitude of the pulsating in- 
stability approaches the classical value obtainable from 
the DNSLE with a weak random noise seeding the insta- 
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bility. but for enhanced quantum fluctuations both the 
ampUtude and the amphtude fluctuations of the pulsa- 
tions become significantly larger than in the classical sys- 
tem. 

In Sec. |TT]we formulate the theoretical model, mainly 
the Gross-Pitaevskii equation (GPE) [23| and its discrete 
variant, the DNLSE We use linear stability anal- 

ysis to find the region of interaction strengths and flow 
quasimomenta where the system develops an instabil- 
ity. In Sec. mil we investigate the time evolution of the 
DNLSE, i.e., classical mean-field theory. Although the 
system initially develops an instability, the subsequent 
time evolution shows nearly periodic recurrence to the 
initial state. In Sec. IIVI we review the well understood 
double well system and argue that the dynamical behav- 
ior of the multi-site system is analogous to the two-site 
system, at least in the limit of weak nonlinearity. In 
Sec.|V]wc study the dynamics beyond the classical mean 
field theory using the TWA. We compare various physical 
properties such as the number fluctuations and the over- 
laps of the state of the system in single realizations with 
an ensemble averages in a few quantitative examples. 
The basic phenomenon of pulsations survive quantum 
corrections, but substantial modifications to the classical 
picture emerge when the atom number is reduced while 
keeping the chemical potential fixed. 

II. THEORETICAL MODEL: DNLSE 

At zero temperature the dynamics of a weakly inter- 
acting BEG in an optical lattice can be modeled by the 
mean field Gross-Pitaevskii equation 0, 

= {'^^ + ^^^^ + ^^'"^'O 

where ^'(x, t) is a wave function corresponding to the 
bosonic field operator such that |\E'(x)p = n(x) equals 
the atom density. The coupling constant 93 is related to 
the scattering length through 53 = Anfi^a/m, where a 
and m are the s-wave scattering length and atomic mass. 
The positive and negative scattering lengths respectively 
correspond to the repulsive and attractive atom-atom in- 
teractions. The Eq. ^ is an approximate description of 
an assembly of a large number of bosonic atoms that are 
in the same quantum mechanical state. Since the equa- 
tion is nonlinear, the coefficient 173 depends on the nor- 
malization of the wave function. Here the wave function 
is normalized to atom number N, so that we have 

J d^x\-^i^)f ^ N . (2) 

The GPE is an increasingly accurate approximation to 
the underlying quantum field theory, for instance, in the 
formal limit ^ 00, 173 ^ with Ng^, and the system 
volume held constant. 



We consider the external potential, V(x), of the form 

l/(x) - ^miivlx^ + Lolrl) + Vo sh? (^^^ . (3) 

Here Vq and dL{= A/2) are the depth and the periodicity 
of the optical lattice. If the harmonic confinement is 
much stronger in the transverse than in the longitudinal 
direction (ujj_ ^ Ux) the GPE can be transformed into a 
one-dimensional form 

with an effective atom- atom interactions gi ~ 2ahLu±. 
Furthermore, when the depth of the optical lattice is 
much larger than the chemical potential of the atoms, one 
can employ the tight-binding approximation. By express- 
ing the condensate wave function ip{x) as a superposition 
of Wannier functions localized within each potential well 
of the lattice, one may obtain the tight-binding version 
of the GPE known as the discrete nonlinear Schrodinger 
equation (DNLSE) 0], 

z^V'/ = -J{i'i+i + ^i-i) + {Vi+x (5) 

We employ periodic boundary conditions. This means 
that for L lattice sites numbered asZ=:0,l,...L— 1, the 
sites I = L actually refers to the site I — and I = —1 
to the site I = L. Physically this corresponds to a ring 
lattice. For other types of boundary conditions there 
obviously have to be some changes in the results. For 
instance, in case of hard- wall boundary conditions a state 
of the atoms that would propagates around a ring will 
instead refiect back from the ends of the lattice. However, 
such variations of the theme will not be discussed further 
in this paper. 

The parameters J and Vi characterize the tunneling 
rate of the atoms from site to site and the external trap- 
ping potential, respectively, whereas x is proportional to 
the atom-atom interactions. Henceforth we ignore the 
non-lattice potential in the direction of the lattice, and 
correspondingly set V; = 0. The value of the interaction 
parameter x depends on the normalization of the wave 
function, i.e., of the complex numbers ipi^ 

^-El'/''!'- (6) 

Differently from both the GPE ^ and Ref. here 
we normalize the state to the one, Af — 1. Even though 
the numbers li/j/p then equal the fractions of the atoms 
residing at the sites I, in what follows wc nonetheless call 
them populations, or downright atoms numbers. The in- 
teraction coefficient is expressed in terms of the (actually 
three-dimensional) unit-normalized wave function of an 
atom in one (actually three-dimensional) potential well 
of the lattice w(x) as 

Anha /" ,3 , , m4 N 
g ^ ^ J d^x\w{^)\\ X=j^9, (7) 
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which expHcitly shows the difference between the inher- 
ently atomic coupling constant, g, and the coupling con- 
stant appearing in the DNLSE, x- Wc frequently re- 
gard the atomic coupling g as variable, which might 
be achieved by employing a Feshbach resonance. The 
DNLSE becomes increasingly accurate, for instance, in 
the asymptotic limit N oo with x L held con- 
stants. Unless it is explicitly stated otherwise, we as- 
sume here that J > and the atom-atom interaction is 
repulsive, x ^ 0. 

In the absence of the remnant trapping potential and 
nonlinearity Eq. ([5]) may be solved with a plane-wave 
ansatz 



J[pl-uj{p)t] 



giving the dispersion relation 



oj{p) = — 2Jcosp . 



(8) 



(9) 



The boundary conditions quantize the lattice momentum 
p to the values 



P 



27r 

T 



-p 



(10) 



where P is an integer that may be chosen to lie in the 
interval [—■§■, ■§ ) • For notational convenience we always 
take the number of lattice sites L to be even. 

When the interaction is switched on, the constant- 
amplitude plane waves are still stationary solutions to 
Eq. IS]) but with a modified dispersion relation 



X 

uj{p) = — 2Jcosp+ — . 

1j 



(11) 



Besides these extended- wave solutions, the nonlinear sys- 
tem may also admit stationary solutions, of the form 
ipi{t) = V'i(0)e~*''*, which are locahzed in space [1^. 
These solutions, so-called gap solitons, usually have an 
energy that lies outside of the linear band spectrum. In 
the continuum model solitonic solutions of the nonlin- 
ear Schrodinger equation can be found in closed form by 
using the inverse scattering method [3(| . However, the 
discrete system has fewer constants of the motion and is 
not integrable as such, so that one has to rely on numer- 
ical techniques for exact solutions. 



A. Modulational instability 

To study the stability of a the plane wave solution ([5]) 
and ([TT]) . we introduce a perturbation around the steady 
state [113, 

^i{t) = ^0(t)[l + u,e'(''-"'*) + t-Je-'^?'-"^*)], (12) 

where q and Clq are the quasimomcntum and the fre- 
quency of the small excitation relative to the initial un- 
perturbed steady state, and Uq, Vq are assumedly small 



mode amplitudes. By the periodic boundary conditions 
the possible excitation modes q are also quantized and 
may be characterized by an integer Q E [— L/2, L/2), so 
that 



q 



2ttQ 



(13) 



(5 = 0, though, corresponds to a change in the normal- 
ization of the unperturbed state ip^ ^ and does not qualify 
as an excitation. 

After inserting Eq. in Eq. ^ with = 0, and 
expanding to the lowest nontrivial order in Uq and Vq we 
get the following matrix equation. 



(14) 



where Uq is the vector [u^, Vq\'^ and is a 2 x 2 matrix 
with the elements 

Mil 



- +4J sm- sm(- +p), 



AI22 = --^ - 4Jsm- sm(--p). 



-Mil 



X 
L 



(15) 



The eigenvalues of M give the small-excitation frequen- 
cies, 



flq ~ 2Jsinpsinq ± Fg; 



J cos p sin' 



2 Q 



2Jcos(p) sin^ I + y 
Z 1j 



(16) 
(17) 



whereas the eigenvectors give the corresponding mode 
amplitudes. 



\ 



Z 1j_ 

2F„ 



\ 



4j cosp sm 7: T i- q + Y 
2F„ 



(18) 



(19) 



As long as the quantity Tq, Eq. p7|l is real, the fre- 
quency of the excitation mode fig is also real and the 
amplitudes Uq and Vq are chosen to satisfy the normal- 
ization conditions \uq\'^ — \vq\'^ = ±1, the sign depending 
on the mode. The reason for this normalization is that 
in quantum theory a unit-normalized state 



Mp, q) = ^(u,e^(f+'')' + z;;e^(f-')') (20) 

then represents one bosonic elementary excitation of the 
condensate. There are two excitation modes correspond- 
ing to the ± signs in Eqs. P^ - P^ for each quasimo- 
mcntum but only those with positive normalization 
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FIG. 1: Coefficient of exponential gain Gq for an excitation 
mode as a function of the parameter Q related to quasimo- 
mentum q hy Q = Lq/{2n). The parameters characterizing 
the interaction strength are A — 0.32 (a) and A = 5.0 (b). 
Modulational instability is only possible for Gq > 0, and ex- 
citation modes may only occur at nonzero integer values Q. 



give physical small-excitation modes. The modes corre- 
sponding to negative normalization are artifacts of the 
present Bogoliubov type analysis, and are henceforth ig- 
nored. For more details on the linear stability analysis 
and its relation to quantum field theory consult, e.g., 

Rcfs. [HI,!!!. 

Since we are dealing with the repulsive atom-atom in- 
teractions, X ^ Oi ^-ll small-excitation frequencies flq 
are real for \p\ < tt/2. However, the existence of com- 
plex eigenvalues cannot be ruled out in the interval 
\p\ S (^, tt], depending on the values of J, x, q and p. If 
for a given p there exists a mode q such that the quantity 
Tq is imaginary, there is a small excitation that grows ex- 
ponentially. This signals a dynamical instability of the 
steady-flow mode p. For a dynamically unstable excita- 



tion mode q we have \uq\ — \vq\ = 
modes cannot be normalized to one 
Figure [T] shows the gain coefficient. 



0, so that unstable 



Ga 



2J 



for two cases (a) A = 0.32, and (b) A = 5.0, where the 
dimensionless interaction strength A is defined to be 



A 



2J 



(22) 



Here we consider a lattice with L ~ 32 sites for p = tt. 
Figure [ija) reveals a single pair of unstable mode with 
Q = ±1, whereas Fig. [Hb) shows four pairs of unstable 
modes. 

Analytically, for any \p\ greater than tt/2 and in the 
limit of large L, the eigenfrequency for the excitation 
mode Q will be complex if 



A > I cos(p)| 



tt'Q' 



(23) 



The corresponding expression for the gain coefficient is 



_ 27tQ^\ cosp|(£A - TT^Q^I cos^ 



and the characteristic time scale for the dynamical insta- 
bility Tq is given by 



2jTa 



1 

Gn 



(25) 



(21) 



The flow with the quasimomentum p is dynamically 
unstable if the inequality ([23|) is true at least for the 
longest-wavelength excitation with Q = I. The critical 
interaction strength approaches zero when the number of 
lattice sites L tends to inflnity, whether for a flxcd atom 
number or if the atom density cx N/L is held constant. 
In this way any flow with \p\ > tt/2 will turn unstable 
with L — > CX). The imaginary part of the complex eigen- 
frequency as well as the corresponding eigenvectors are 
the same for the modes Q and —Q [see Eqs. ([T8| . ([T9|) ]. 
so that we regard these two modes as equivalent as it 
comes to the instability. 



III. TIME EVOLUTION AND PULSATING 
INSTABILITY 

We carry out numerical simulations on the DNLSE to 
study the growth of the unstable excitation mode in a lat- 
tice for a suitable range of interaction parameters. For a 
given number of lattice sites L and flow momentum p, the 
number of unstable modes in the linear stability analysis 
depends only on the interaction parameter A, Eq. (j22p . 
Here we focus on long- wavelength excitations in the limit 
of weak atom-atom interactions. Two numerical methods 
have been used for the time evolution, an unconditionally 
stable Crank-Nicholson type algorithm [33| and a sixth- 
order accurate FFT split operator algorithm that works 
in the same way for DNLSE as is discussed in Rcf. [s^l 
with the ordinary GPE. 
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FIG. 2; (Color online) Collapse and revival of the density 
pulse in the BEC evolution when the lattice starts from a 
dynamically unstable state in a density plot for the popula- 
tions of the sites the lighter shades corresponding to 
larger populations. The parameters are A = 0.48, p = tt, 
and L = 32. These parameters correspond to a single un- 
stable mode, as per linear stability analysis. The instability 
drives the initial homogeneous atom distributions into a den- 
sity peak that subsequently disperses back to the initial state, 
and the process repeats. 



A. Single unstable mode 

A perusal of the condition shows that the range of 
interaction strengths where the Q — \ mode is unstable 
but the Q — 2 mode is not is given by 



2 

I cospl — < A < 4| cosp 
ij 



L 



(26) 



Figure [5] shows a typical density plot of the time evo- 
lution of the BEC initially prepared in the plane wave 
state at the edge of the Brillouin zone {p = tt) seeded 
with random Gaussian noise, for the number of lattice 
sites L — 32 and the interaction strength A ~ 0.48. This 
value of A corresponds to one unstable mode in the lin- 
ear stability analysis. Although a tiniest amount of noise 
(either in real experiments or in numerical simulations) 
is sufficient to trigger the instability, an external noise 
of amplitude ^ 10~* is added to speed it up. The same 
applies to all of our demonstrations of the pulsating insta- 
bility in multisite lattices up to Sec.|Vl where we adopt a 
model for the quantum noise instead. It has been tested 
in a number of runs that the time for the onset of the in- 
stability for fixed values of the other parameters depends 
logarithmically on the amplitude of the added noise, as 
it should. 

Figure [3] depicts snapshots of a pulse that moves dur- 
ing its formation from the initial flow state with quasi- 
momentum p = 15 tt/16. In this figure we take a larger 
lattice with L = 128 sites and the interaction strength is 
A = 0.25, so that there is still one and only one unstable 
mode. The pulsating behavior of the peak still persists 
but the peak moves. By virtue of the periodic bound- 
ary conditions a pulse that goes over the right edge will 
reappear at the left edge of the lattice. 
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FIG. 3; (Color online) Density evolution of the BEC for the 
initial flow momentum different from tt. Here the parameters 
are L = 128, A = 0.25 and p = 157r/16. 
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FIG. 4: Fraction of the initial state f{t) in the state of the 
lattice plotted as a function of time. Here the parameters are 
L = 32, A = 0.48 and p = tt. 



The pulsating behavior also manifests in the fraction 
of the initial state ipiiO) remaining in the state of the 
lattice at a given time. 



m 



(27) 



In Fig. m we plot the overlap, f{t), as a function of time 
for the same parameters as in Fig. [2l It is revealed that 
the instability drives the system far from, and subse- 
quently brings it back to, the original unstable steady 
state, and the process repeats. Each dip in the plot rep- 
resents formation of a pulse during the course of time. 
It is also noted that the quantity f{t) does not vanish 
all the way to zero, indicating that the pulsed state is 
not orthogonal to the initial steady state. Furthermore, 
a closer inspection of this plot shows that the subsequent 
peaking events are not strictly periodic; the interval be- 
tween the dips varies slightly, implying a quasi-periodic 
phenomenon. 



By analyzing data of this kind a number of obser- 
vations emerge. First, in contradiction to the common 
behef that the instabihty may develop into an irregu- 
lar dynamics, here it causes the atoms to pile up in a 
single-peaked distribution of the site populations 
Furthermore, upon continued time evolution, the system 
returns very close to the initial unstable state, again pul- 
sates to a peak, and so on. We have nearly periodic peak- 
ing and recurrences to the unstable initial state. Second, 
for the initial flow state p ^ tt, the pulse also moves 
with a velocity that turns out to be close to the group 
velocity of the carrier wave, Vg = 2Jsinp. Third, the 
peak may occur at any lattice site. It is the random 
noise that seeds the position of the peak. This is in ac- 
cordance with translational invariance of the lattice: A 
lattice-translated pulsed solution is also a solution of the 
DNLSE, and there is no preferred lattice site for the oc- 
currence of the pulse. Fourth, for weak initial noise the 
dependence of the pulsation phenomena on the initial 
noise is otherwise weak. For instance, even if the timing 
of the first pulse depends logarithmically on the noise am- 
plitude and the spacing of the subsequent peaks is also 
affected, the dominant time scale is evidently the time 
scale of the instability from Eq. (|25p . In Fig.[4]the peak 
spacing equals approximately 10 r,. 
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FIG. 5: (Color online) Density plot of the populations in 
the case where there are four unstable modes. The parameters 
are L — 128, A = 2.0, and p — tv. Lighter shading represents 
higher site populations. 



B. Multiple unstable modes 

For a lattice with a large number of sites the one-peak 
condition (|26p is highly impractical because the interac- 
tion strength A needs to be extremely small and the pulse 
revival period is correspondingly long. For instance, the 
shortest possible time scale for the dynamical instabil- 
ity in the one-peak case, r, is determined from Eqs. (j26p 
and 1^ as 2Jr = L"^ / {2V3tt^ \ cos p\) oc L^. 

On the other hand, the condition that exactly Q modes 
are unstable is obtained from Eq. (|26p in the form 

2 2 

g2|cosp| ^<A<(Q + l)2|cosp| ^. (28) 

For reasonable interaction strengths there might be more 
than one unstable mode. Figure[5]is a typical representa- 
tive of the dynamics of the BEC with multiple unstable 
modes. Here we take A'' = 128, A = 2.0, and p = tt. 
Each bright spot represents a pulse. As before, the right 
edge of the plot wraps around to the left edge by virtue 
of the periodic boundary conditions. Equation (pS)) gives 
five unstable modes, while most of the time we see three 
or four pulses in this example. However, these pulses 
are not independent. Presumably because of nonlinear 
mode-mode interaction, they move around, join and split 
as they collapse and revive. 



C. Evolution in Fourier space 

The recurrences observed in the peaking events in the 
DNLSE resemble the energy recurrences in the Fermi- 
Pasta-Ulam (FPU) problem [H. The FPU model deals 
with the evolution of a lattice chain with nonlinear in- 
teractions between the nearest-neighbor atoms when ini- 
tially a single low-energy mode is excited. For a time 
scale much longer than the time period of the normal 
modes, the energy is well localized to the given excited 
mode, while the amplitudes of the higher-energy modes 
decay exponentially as a function of the energy difference 
from the initially excited mode. For a longer time scale it 
has also been noticed that recurrence of the initial excited 
mode is possible. 

The pulsating behavior of the density distribution of 
the BEC atoms in the lattice can be viewed as a similar 
recurrence phenomenon as observed in the FPU model. 
We have started with a steady state for a given flow quasi- 
momentum (p > tt/2) and a suitable nonlinear interac- 
tion strength to trigger the instability in the system. Er- 
godicity immediately suggests that the energy initially 
fed into a single mode should distribute evenly between 
all Fourier modes. However, the excitation amplitudes 
of the modes other than the mode corresponding to the 
initial steady state seems to decay exponentially with the 
index Q. The energy localization to a few Fourier modes 
in a nonlinear system has been studied also in other sys- 
tems [1^. The existence of discrete breathers in a non- 
linear lattice system is an example. Recently, energy lo- 
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FIG. 6: (Color online) Time evolution of the Fourier modes 
of the lattice. The parameters are L = 32, p = tt and 
A = 0.48. The mode Q = corresponds to the initial steady 
state while modes Q = 1,2,3 .. . are the low-lying excitations. 
Only the modes with Q > are shown. 



calization in Fourier space in a so called 'q-brcather' has 
been investigated in [36j . 

Figure [5] shows the time evolution of the Fourier modes 
-0, of DNLSE for the parameters A = 0.48, = 32. 
Only a few components are seen to be excited, as the 
amplitudes of the higher-energy modes are suppressed 
exponentially. The dominant Fourier components are P 
and P ± 1 which implies that the excitation modes that 
go unstable should have the indices Q = ±1 with respect 
to the initial steady state, as expected. 



IV. DOUBLE WELL ANALOGY 

In order to explain qualitatively the pulsating behav- 
ior of the density distribution of the BEG in the lat- 
tice, we study a coupled double-well system. Writing 
"00,1 = 100,1 1^*'^"'^ the dynamics of such a system can be 
described by a pair of equations (37j 



z{t) = -VT^2^sin<?!)(i) 
z{t) 



(j){t) = Az{t) + 



cos (f>{t). 



(29) 



where z = IV'ip — IV'oP and (j) = — are the fractional 
population imbalance and the relative phase between the 
two wells. The normalization again is l^oP + IV'iP = 1- 
The correct evolution for z and (j) is obtained from the 
Hamiltonian 



H 



Az2 



(30) 



by regarding these quantities as canonical conjugates. 
The Hamiltonian (energy) is a constant of the motion 
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FIG. 7: Representative constant-energy contours of the two- 
site Hamiltonian (|3Up in the {cj>, z) plane for the interaction 
strengths A = 0.5 (a) and 1.5 (b). 



and the double-well system is therefore integrable. The 
evolution of the system keeps it on a constant-energy 
contour in the space with z and (f) as the axes. 

By inspection it can be checked that the fixed points 
of Eq. ((29|) are 2 = 0, = mi, where n is an integer. A 
potentially unstable steady state in the multiwell system 
corresponds in the two-well system to the steady state 
z ~ Q, (j) ~ TT. The behavior of the orbits near this equi- 
librium point can be examined by using linear stability 
analysis as before. It can be easily verified that the state 
(0,7r) is stable for the values A < 1, and unstable other- 
wise. 

In Fig. [7] we show constant-energy contours of the two- 
well system for A = 0.5 (a) and 1.5 (b). We have drawn 
the (p axis from to 27r so that the potentially unstable 
fixed point (0, tt) lies at the centers of the panels. For 
A = 0.5 the potentially unstable steady state is an ellip- 
tic fixed point and the time evolution takes the system 
periodically around this point. At A = 1 the elliptic fixed 
point bifurcates, and for A = 1.5 there is a homoclinic 
orbit with the emergence of two symmetric off-centered 
elliptic fixed points. Thus, starting in the vicinity of 
what used to be the potentially unstable steady state, 
the system takes off in an unstable direction along the 
homoclinic orbit and goes around one of the bifurcated 
elliptic fixed points. 

The double-well system also admits an analytic solu- 
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FIG. 8: Time evolution of the population imbalance z given 
by an exact analytical solution of Eqs. (|29|) for the parameter 
A = 1.5 and the value of the Hamiltonian H = 1.0. 



tion in terms of Jacobian elliptic functions. In Fig. [5] we 
plot such a solution for z[£) for the parameter A = 1.5 
and the conserved value of the Hamiltonian H — 1.0. 
The double-well system is then unstable in linear sta- 
bility analysis and the energy corresponds to an energy 
contour close to the homoclinic orbit, see Fig.lHt^b). The 
oscillations in the population imbalance bear a striking 
resemblance to Fig. 2] with the plot of the overlap in the 
multiwell lattice, and may be viewed as an analogue of 
the pulsating instability. 

We have presented two different views on the pulsating 
instability in the two-site system, energy contours and 
time evolution of the population imbalance, but they de- 
scribe the same physics. The dynamically unstable sys- 
tem performs periodic oscillation where the system re- 
peatedly recedes far away from the unstable state and 
subsequently returns. 

The multisite system basically shares the behavior of 
the two-site system in a multi-dimensional phase space: 
Starting from random noise in the neighborhood of an un- 
stable steady state, it evolves away from, and returns to, 
the initial state, and the process then repeats. These pe- 
riodic recurrences occur in a 2L-dimensional phase space 
on a constant-energy surface in full analogy with the two- 
site system. The two-site system is strictly periodic since 
there is no motion out of the constant-energy surface, a 
one-dimensional curve. However, in the multiwell case 
the dimension of the constant-energy surface is 2L — 1 . 
Our pulsating instability strongly suggests that the mul- 
tiwell system stays close to the homoclinic orbit while it 
evolves. However, depending on the initial noise it still 
has a large state space to explore, and the motion may 
deviate slightly from the homoclinic orbit. Upon looping 
around one of the stable fixed points, the multisite sys- 
tem therefore does not have to return to exactly where it 
started from. This might account for the slight variations 
in the period of the pulsations. 



In short, we attribute the pulsating dynamics in the 
multiwell system to a remnant of integrability. The mul- 
tiwell system evolves along a similar homoclinic orbit 
that exists in a two-well system, and the motion away 
from the homoclinic orbit remains bounded. 



V. TRUNCATED WIGNER APPROXIMATION 

The GPE can constitute a very accurate modeling of a 
weakly interacting BEG. In Sees. HIHTVl wc discussed the 
dynamics of a BEG within the classical mean-field the- 
ory, including double- well and multi-well systems. We 
know in the case of double-well systems, however, that 
when the tunneling barrier becomes large, quantum and 
thermal fluctuations can considerably influence the atom 
dynamics (ssj . Similarly, in optical lattices the kinetic 
energy is represented by the hopping of the atoms be- 
tween adjacent lattice sites, and the hopping amplitude 
can be significantly reduced in deep lattices. Quantum 
correlations and fluctuations are a matter of balance be- 
tween atom-atom interactions and hopping, and may be 
correspondingly enhanced. In the following we include 
quantum fluctuations in the dynamics in a lattice using 
a stochastic phase space method, the truncated Wigner 
approximation (TWA). 

TWA was introduced for multi-mode dynamics in non- 
linear optics [121. Details how to implement TWA 
in different atomic BEG systems may be found, e.g., 
in Rcfs. [H, m, m, m, studies of zero- and finite- 
temperature nonequilibrium dynamics in ID lattice sys- 
tems are documented in Refs [l5|, [l^, [1^, [H, HO, HTj , and 
the effects of dynamical instabilities in lattices have been 
exphcitly addressed in Refs. [H, [H, [3|. In the TWA 
one neglects the third-order derivatives in the general- 
ized Fokker-Planck like equation for a Wigner distribu- 
tion function [43| . This permits a nonlinear stochastic 
differential equation for a classical field that represents 
the quantum field. Here the equation is just the DNLSE. 

In our present case of TWA quantum fluctuations and 
correlations enter only through the stochastic initial state 
of the classical field. Our emphasis on quantum fluctu- 
ations is different from typical finite-temperature dom- 
inated TWA approaches in higher dimensions [13] ■ In 
an attempt to capture quantum mechanics as accurately 
as practicable, we therefore pick the initial states of the 
wave function using the Bogoliubov approximation. 

In our numerical simulations wc consider a stable sta- 
tionary superfluid flow that is instantaneously, at t = 0, 
rendered dynamically unstable by changing atom-atom 
interactions. We investigate the effect of quantum fluctu- 
ations by considering two different types of initial states: 
a noninteracting BEG with all atoms in the same one- 
particle state, and an interacting system in which the 
atom-atom interactions force fraction of the atoms out 
of the condensate. We first generate a stochastic initial 
state = 0) accordingly, then the TWA dynamics 
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follows from the DNLSE 



with 



. d w 



w 



(31) 



This process is repeated for a number of initial states. 
Two types of results are of interest: Individual trajecto- 
ries ipj^ {t) represent the outcomes of individual exper- 
iments, and appropriate averages over the collection of 
trajectories are used to calculate pertinent quantum ex- 
pectation values. The present TWA formalism is very 
similar to the one used in Rcfs. [IE, H^, E^] , except that 
in each TWA realization we fix the total atom number 



A. Initial state 

In order to describe the generation of the stochastic 
initial states for the TWA dynamics we begin with the 
quantum version of the Bogoliubov approximation. We 
again consider the stationary solution for a moving plane 
wave as in Eqs. ^ and ([TT|) . The fluctuations around the 
stationary solution are governed by the decomposition of 
the atom field operator 

Mt) = ^°{t)(^o + SM-t), (32) 

where the total number of condensate atoms, 

is assumed much larger than one and Sipi {t) is supposedly 
"small." Heretofore we make a difference between total 
number of atoms iV, number of condensate atoms Nc, 
and number of noncondensate atoms iV„. Analogously 
to our earlier classical Bogoliubov treatment, the fluc- 
tuation part of the atom field operator, 5il}i{t), can be 
written in terms of quasiparticle operators ctq and a^^ as 

(5^i(t) = ^^(uga,e['(f+«)'-'^^'*l+w;dte*[(f-9)'+*f^^*l). 

(33) 

Within the Bogoliubov approximation the quasiparticles 
make a bosonic ideal gas, or alternatively, the operators 
Ciq may be viewed as the lowering operators for a col- 
lection of independent harmonic oscillators. The normal 
mode frequency f2g and the quasiparticle amplitudes Ug, 
and Vq are given by Eqs. (fTH]) . (fTS]) . and (fT9|) : only the 
excitation modes with positive normalization arc used. 
The unit normalizations for both the unperturbed plane 
wave ij/^ and the amplitudes of the excitation modes Ug, 
Vq were fixed earlier so that these quantities work as given 
also in the present quantum version of Bogoliubov theory. 

The expectation value of the number of non- 
condensate atoms in the Bogoliubov theory is 



E(KI' 



(34) 



{a\aq) 



At T = we have (atdg) = 0, and the non-condensate 
atom number is simply 



(36) 



In order to construct the initial state of the lattice 
for TWA evolution we replace the quantum operators 
(dg,dj,) in Eq. ^6'6\ by complex stochastic variables 
{uq, a* I ) obtained by sampling the corresponding Wigner 
distribution function. Our formalism follows Ref. \2&i. 



except that here we fix the total atom number [4J|. The 
Wigner function at T = is [isj 



Wiaq 



■exp[-2|ag|2] 



(37) 



The function W{aq,a*) is a Gaussian and we easily 
find the average number {a*aq)w = 5 of excitations 
in each mode q owing to quantum noise. Each unoc- 
cupied excitation mode exhibits uncorrelated Gaussian 
noise, distributed over the plane wave basis, and normal- 
ized to an average of a half an excitation quantum per 
mode. The attendant atoms provide seeding for scat- 
tering events in the TWA dynamics. However, Wigner 
functions correspond to symmetrically ordered expecta- 
tion values of functions of the operators ctq and dj, and 
in the end the expectations values {a*aq)w = ^ need to 
be subtracted to obtain normally-ordered quantum av- 
erages. For each stochastic realization, the number of 
non-condcnsatc atoms therefore reads 



(38) 



This may fluctuate about the mean value 7V„ = J2q bgP 
for each realization, but correctly averages to Nn- 

Since the total atom number TV is fixed, the number of 
condensate atoms in each individual run is given by 



(39) 



which also fluctuates around the mean value Nc ^ N — 
Nn- In the initial state for each realization we then set 



ao = ^Nc + 1/2 , 



(40) 



so that, in the place of Eq. ([32|) . wc have the initial state 
for the TWA method 

(0) - ^'?(0)ao+4= ' 

(41) 

where aq and a* are stochastic variables in each realiza- 
tion of the TWA. Note that, even though we consider 
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a uniform system with a plane wave phonon basis and 
uncorrelated noise in the initial phonon modes, the fix- 
ing of the total atom number introduces long wavelength 
correlations in the system between the condensate mode 
and the excited quasiparticle modes [iij . 

Based on the Bogoliubov results it is interesting to es- 
timate when the classical description of the condensate 
dynamics from the DNLSE. as in Eq. approaches 
the quantum results, as per the TWA. For instance, if 
the ratio of non-condensate particles to condensate par- 
ticles in the initial state is large, one may expect the 
classical DNLSE to provide an especially poor approxi- 
mation. Now, the expression (jl9p for the amplitudes Vq 
shows that the interaction strength enters through the 
dimensionless parameter %/ J. By virtue of our normal- 
izations, the interaction parameter x, Eq. ([7]), is propor- 
tional to the number of atoms. If the number of atoms is 
increased while adjusting the strength of the atom-atom 
interactions in such a way that x/ J is held constant, the 
number of noncondensate atoms remains constant and 
the ratio of non-condensate particles to condensate par- 
ticles in the initial state is reduced. We expect the error 
in DNLSE correspondingly become smaller [l^. In fact, 
even if the non-condensate population in the initial state 
vanishes, the relative contribution of the vacuum Wigner 
noise is reduced when the condensate population is in- 
creased, and we again expect the error in DNLSE to be 
reduced. 

The potential instability, however, adds a new element 
to the picture. Namely, for the consistency of our present 
development the excitation modes must be normalized to 
— = 1, whereas for an unstable mode necessar- 
ily \uq\'^ — = holds true. This means that when 
the point of dynamical instability is approached continu- 
ously, for some mode q the amplitudes diverge, \uq\ — *■ oo 
and \vq\ — > oo. By Eq. (|36p the number of noncondensate 
atoms also diverges, Nn — *■ oo. At the point of instability 
steady-state quantum fluctuation diverge, according to 
the Bogoliubov approximation. The Bogoliubov method 
is based on the assumption of small fluctuations so that 
it in itself fails for an unstable system, but it is clear that 
the onset of instability is associated with large quantum 
fluctuations. 

While the asymptotic results are easy to infer, only 
a detailed analysis can tell what the number of noncon- 
densate atoms is in a specific situation. Figure [9] shows 
the number of noncondensate atoms Nn as a function 
of the interaction parameter A for the number of lattice 
sites L = 32. Both the generally very small number of 
noncondensate atoms for weak interactions and the di- 
vergence at the onset of instability A = 0.308 ~ tt^/L 
are evident. 



B. Numerical realization 

We study the non-equilibrium quantum dynamics of a 
BEC within TWA. We consider a BEC in a lattice that 
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FIG. 9: Number of noncondensate atoms iV„ given by Eq. p6|) 
as a function of the scaled interaction strength A within Bo- 
goliubov approximation for L — 32. Note that the critical 
interaction strength at the onset of the instability is 0.308. 



is initially in a stable steady flow state. This is crucial 
for the validity of the Bogoliubov approximation for the 
initial state as well as for the generation of the stochastic 
initial noise in the TWA. A stable state of the atomic gas 
can be experimentally driven to a dynamically unstable 
state, for instance, by starting with a moving lattice with 
the quasimomentum of the atoms p < tt/2 and accelerat- 
ing the lattice so that the momentum satisfies p > tt/2. 
In our examples, though, we envisage suddenly modifying 
atom-atom interactions. Henceforth we always have the 
number of lattice sites L = 32 and the initial flow state 
with p ~ TT. We perform the simulations for two differ- 
ent types of initial states: noninteracting atoms with all 
atoms in the condensate, and an interacting state with a 
nonzero condensate depletion. As to the interacting ini- 
tial state, we select the dimensionless interaction strength 
A = x/2</ = 0.284. For this value the average noncon- 
densate atom number is iV„ ~ 30 whereas the critical 
value of A at the onset of the instability is 0.308 , as per 
Fig. O In all cases considered here, the initial depletion 
Nn/N is less than 10%, and the Bogoliubov approxima- 
tion for the initial state is presumed accurate. At the 
beginning of the time evolution the atom-atom interac- 
tions arc turned up instantaneously to the dynamically 
unstable value of A = 0.48. We investigate the effect of 
quantum fluctuations by varying the total atom number 

and the atomic interaction strength g (Eq. ([71)) so that 
the interaction parameters x and A (and the chemical po- 
tential) remain constant. 

For each individual realization of the time evolution of 
the ensemble of the Wigner distributed wave functions 
we sample the initial state as explained in the previous 
Section [VB I The random numbers (ctg, a*) are obtained 
using the Box-Mucller algorithm (3^. As before, we inte- 
grate the dynamical equations PT|) using the FFT split- 
step method [3^ . Depending on the task at hand, we 
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either study the individual trajectories ipj^{t), or aver- 
ages over a number of trajectories. 

It should be noted that our process contains an uncon- 
trolled approximation. There are some indications [isl . 
[26l . [45I . [ia i that time evolution as prescribed by TWA 
could becomes exact in the asymptotic limit N 00 
with a fixed x- but even if this were the case, the asymp- 
totic limit per sc does not tell how good the results are for 
given values of the parameters. For instance, we do know 
how long the TWA time integration remains valid. This 
question should be studied fully quantum mechanically, 
but at the present time we do not know of any framework 
to address the issue. How would one describe the modula- 
tional instability of a classical nonlinear system ab initio 
using linear quantum mechanics, which by definition pre- 
scribes a stable evolution of any initial state? Moreover, 
in a classical instability the translation invariancc of the 
state of the lattice is broken, but unitary time evolution 
from quantum mechanics will not spontaneously break 
any symmetry of the Hamiltonian. For the time being 
we rely on the practical observation from working with 
the TWA method that it usually signals its own demise. 
Typically, upon integrating long enough, each individ- 
ual trajectory seemingly loses any relation to the initial 
state. In our computations we have not encountered a 
situation of this kind. 



£ 0.5 



£ 0.5 





i 






1 




- 

(a) 















1000 2000 
2Jt 


3000 


\ 




\ 


\ \ 


( 




(b) 






1 







S 0.5 - 




2Jt 




1000 2000 3000 1000 2000 3000 
2Jt 2Jt 



FIG. 10: (Color online) Overlap of a single stochastic real- 
ization with the initial state as a function of time sampled 
from the Wigner distribution for a fixed chemical potential 
with different atom numbers (a) = 10*^, (b) iV = lO", (c) 
N = 10^ and (d) A'' = 500. The initial state is interacting with 
Nn = 30. The parameters of the simulation are A = 0.48, 
L = 32, and p = n. There is no significant damping in the 
oscillation even if the noncondensate noise is substantial. 



C. Results 

Since the Wigner functions govern symmetrically or- 
dered expectation values, we need certain transforma- 
tions to calculate the ordinary normally ordered expec- 
tation values from the simulation data |26l . l40j . Here we 
only consider the lowest-energy band in the tight-binding 
approximation, so normally ordering the operator expec- 
tation values is straightforward. According to Ref. [2q . 
we have the atom number in a lattice site 

ni = {i:^4,i)w-\, (42) 
with the corresponding fluctuations 



An; = ^j{{rii'i?)w - {rii'i?w ~ 4 • (43) 

The overlap of the field amplitudes ^/'/(O) and ^/';(i), which 
is a measure of the revival of the pulse, is given by 

r{t)^{\i:iri{m0f)w. m 

In Fig. [To] we show typical single-trajectory results for 
the overlaps f{t) for different values of the atom-atom 
interactions g and the total atom number N . All dif- 
ferent cases have the initial interaction parameter A = 
Ng/{2J) = 0.284 that is instantaneously changed to 
A = 0.48, but the values of N and g are varied individ- 
ually so that the total number of atoms is (a) N = 10^, 
(b) N = 10'^, (c) N = 10^, and (d) N = 500. 



The number of noncondensate atoms in the Bogoliubov 
approximation docs not change whenever we vary Nc and 
g for a fixed value of Ncg. For a smaller total number of 
atoms N, the noncondensate atom fraction Nn/N in the 
initial state and the coupling constant g are larger and, 
consequently, quantum effects are stronger. Each plot 
in Fig. [TO] clearly shows the pulsating instability with- 
out any noticeable damping. Quantum fluctuations and 
the initial noncondensate population act as a seed for 
the scattering processes taking atoms out of the ground 
state. The period of the pulsating instability depends on 
the interaction coefficient g, albeit not dramatically; the 
larger the value of g (corresponding to the smaller value 
of N), the shorter the period of oscillation, since higher 
scattering rate leads to a faster ground state depletion 
and, consequently, a shorter oscillation period. 

Figurc[TT]rcprcscnts an ensemble average of the overlap 
f^{t) sampled over 400 trajectories for the total number 
of atoms (a) N = W^, (b) N = 10^, (c) N = 500 and 
(d) N = 300 for the interacting (solid curve) and the 
noninteracting (dashed curve) initial states. In all cases 
the final value of the interaction parameter is again A = 
0.48 (oc Ng), but the nonlinearity g and the total atom 
number are different. 

While the almost undamped quasiperiodic behavior 
is seen in individual stochastic realizations, the quan- 
tum mechanical ensemble averages of the wave function 
revival become progressively weaker when the effective 
interaction strength is increased. This is because the 
shape and the timing of the pulsations in each realization 
change due to quantum fluctuations. Enhanced quan- 
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FIG. 11: (Color online) Comparison of the ensemble average 
of the overlap of the state of the lattice sampled over 400 
realizations for the initial number of non-condensate atoms 
iV„ = (solid line) and iV„ — 30 (dashed line), and for the 
total number of atoms (a) iV = lO", (b) iV = 10^, (c) = 500 
and (d) A'^ — 300. The parameters of the simulations are 
A = 0.48, L = 32, and p = n. The two different initial states 
represent similar time evolution when the total atom number 
is large. However, the curves start deviating as the number 
of atoms is reduced. 



turn fluctuations for smaller N (and larger g) are clearly 
observable in increased damping rates. For the interact- 
ing initial state, the initial noncondensate atom number 
Nn = 30 provides a nonnegligible noncondensate atom 
fraction only if the total atom number is small, and so 
the results for the noninteracting and interacting initial 
states differ significantly only for small total atom num- 
bers. Figurc fT2l represents a comparison of the overlaps in 
a typical single realization and in an ensemble average for 
the large atom number N = 10^, and with the interact- 
ing initial state. The other parameters of the simulations 
are same as in Fig. [TTl 

Similar damping is also observed in the atom number 
fluctuations at a flxed site, say, center of the lattice. In 
Fig. [13] we show ensemble averages of the atom number 
fluctuations An/^yn for the total number of the atoms 
(a) N = 10^, (b) = 10^, (c) TV = 500 and (d) N = 300 
and for the interacting initial state. The parameters of 
the simulations are as before. The pulsating instability 
causes strongly super-Poissonian atom number fluctua- 
tions at any flxed lattice site. 

The main focus in this work is to estimate the inher- 
ent quantum effects on the pulsation phenomenon of the 
BEC in an optical lattice. We have found that quantum 
fluctuations have an effect on the collapse and revival of 
the pulse. In the case of a single realization the revival is 
very robust and repeats over a long time. The ensemble 
averages of the overlap of the wave function revival dis- 
play a recognizable damping behavior due to quantum 
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FIG. 12: (Color online) Comparison of the overlaps in a typ- 
ical single realization (solid line) versus an ensemble average 
sampled over 400 realizations (dashed line) with the initial 
interacting state for iV„ = 30 and = 10^. The parameters 
of the simulations are A = 0.48, L = 32 and p = n. 
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FIG. 13: (Color online) Atom number fluctuations in the cen- 
tral lattice site as a function of time sampled over 400 real- 
izations for the total numbers of particles (a) A*' = 10*, (b) 
A' = 10^ (c) N = 500 and (d) N = 300 with Af„ = 30. We 
use the interacting initial state. 



fluctuations. We have also found that the damping rate 
is increased if the initial state is depleted with a nonneg- 
ligible noncondensate fraction. 

The TWA also provides for computing of quantities 
that are affected by quantum fluctuations and could be 
readily measurable, but may be very difficult to extract 
from the full quantum solution of the problem if it were 
on hand. Parameters for the shape of the pulse, obscured 
in quantum mechanics by averaging over the random time 
and position of the pulse, are an example. In Fig. 1141 we 
show the amplitude for the first pulse in the pulsating 
instability as a function of the total number of atoms 
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FIG. 14: (Color online) Ensemble average of the amplitude 
(left) of the first peak and its fluctuations (right) during pul- 
sation as a function of the total number of atoms, A^, for a 
fixed chemical potential for interacting (solid line) and non- 
interacting (dashed line) initial states. The parameters of the 
simulation are same as in Fig. 1111 The amplitude reaches the 
classical value for higher atom numbers. Quantum noise in- 
creases the amplitude and the fluctuations of the amplitude 
for small atom numbers. 

(for a fixed nonlinearity, as before), and its fluctuations. 
These are averaged over 400 realizations, both for the 
interacting (solid curve) and the nonintcracting (dashed 
curve) initial states. The parameters of the simulations 
are the same as in Fig. [11] 

The pulse ampUtude approaches a classical value that 
may be obtained from the GPE in the limit of weak 
quantum fluctuations (here corresponding to large N and 
small g). However, both the amplitude of the atom 
pulse and the amplitude fluctuations grow when the atom 
number becomes small and quantum fluctuations become 
more dominant. Similarly to the case of the overlap re- 
vival, the distinction between the two initial states is ap- 
parent only when the noncondensate fraction is nonneg- 
ligible. 

VI. CONCLUDING REMARKS 

In this paper we have presented new insights into the 
unstable dynamics of the BEC in an optical lattice in the 
limit of weak atom-atom interactions, with quantum fluc- 
tuations included. The common belief is that the flow of 
the dynamically unstable BEC in an optical lattice would 
be erratic, or lead to the formation of stable solitons. 
Here we move a step further and show that, in classi- 
cal mean-field theory, the instability may also trigger a 
quasi-periodic pulsation in the atom density distribution 
if the atom-atom interactions are weak. The requirement 



that linear stability analysis finds a single unstable mode 
gives the scale for the 'weak' nonlinearity and the ensuing 
pulsating phenomena. 

A qualitative argument has been put forward to ex- 
plain the pulsating behavior of the dynamics by com- 
paring the lattice system with the integrable double-well 
system. In the case of two wells the unstable mode leads 
to a non-trivial dynamics in the population imbalance 
such that an infinitesimal noise could produce a large- 
amplitude collective oscillation of the atoms between the 
wells. We surmise that the pulsating instability is a rem- 
nant of the integrability as manifest in the two-well sys- 
tem. 

We incorporate quantum fluctuations using stochastic 
phase-space methods. We use the Bogoliubov approxi- 
mation to generate the initial state for the time evolu- 
tion of the system. A sequence of the stochastic fields 
obtained in this way are then used to calculate expecta- 
tion values of the observables. We also compare the single 
realization results with the ensemble averages. Generally 
speaking, the pulsating behavior survives in the face of 
quantum fiuctuations. It is observed that the quasiperi- 
odic behavior in the time dynamics can still be seen in 
single realizations. However, the quantum averages show 
that the revivals of the pulses tend to get washed out as 
the atom number gets smaller while the chemical poten- 
tial is held constant. 

For experimental realizations, the fiow states p w tt 
near the Brillouin zone boundary can be prepared by 
accelerating the lattice [Ij. Alternatively, by exploiting 
the symmetry of the DNLSE, for every solution ■4'i{t) 
with the given interaction parameter A there is a solu- 
tion (— l)V/*(i) for —A. That means the state for p = tt 
in the repulsive case is equivalent to the state for p — 
in the attractive case. Manipulation of the sign of the 
interactions gives much additional leeway for the exper- 
iments. However, given that the pulsating phenomenon 
only occurs in the limit of weak nonlinearity, the corre- 
sponding time scales for the pulsation can be very long 
and pose severe technical challenges. 
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